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Abstract 

We generalize the Dogterom-Leibler model for microtubule dynamics |DL] to the 
case where the rates of elongation as well as the lifetimes of the elongating and short- 
ening phases are a function of GTP-tubulin concentration. We study also the effect 
of nucleation rate in the form of a damping term which leads to new steady-states. 
For this model, we study existence and stability of steady states satisfying the bound- 
ary conditions at x = 0. Our stability analysis introduces numerical and analytical 
Evans function computations as a new mathematical tool in the study of microtubule 
dynamics. 
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General Theory We analyze some mathematical aspects of the phenomenon of dynamic 
instability of microtubules. Our theoretical model includes a simplified, semi-infinite ge- 
ometry, in which infinitely rigid microtubules grow perpendicularly to a nucleating planar 
surface (See Figure 2). Dogterom and Leibler studied this theoretical model with this 
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semi-infinite geometry, neglecting concentration variations in the process of assembly and 
disassembly of microtubules [DLJ. In this model infinitely rigid microtubules are growing 
perpendicular to a nucleation planar surface (Fig. 2) and randomly switching between the 
growth state (+) and shrinkage state(— ) or polymerization and depolymerization. 

In this model each microtubule switches randomly between the assembly state +, in 
which it grows with the average speed v + proportional to the local monomer density c, and 
the disassembly state — , in which it shrinks with the average speed v~ . The frequencies of 
the transitions between the two states, /7 (of the "catastrophies" from + state to — state) 
and (the "rescues" from — state to + state), determine, together with v + , v~ determine 
the behavior of the model. We summarize now the main results obtained through Evan 
function and time-evolution simulations. The main result is the existence and stability 
of steady state solutions for a convection diffusion equation modeling the above-described 
process. The novelty in this study is the analysis of the general case when the dynamics 
parameter are depending on local concentration. The mathematical tools introduced to 
handle the more complicated resulting equations should be of general use. 

1 Introduction 

Microtubules are natural polymers found inside of living eukaryotic cells. Each polymer 
subunit is an obligate dimer of proteins synthesized from the alpha and beta tubulin genes. 
When concentrated above a threshold level, tubulin dimers bind to each other forming a 
beta lattice of subunits that bind head to tail and side-to-side. Under appropriate condi- 
tions, the inherent curvature of the lattice sheet produces a tube, typically 13 subunits in 
cross-section, which can extend at the tube ends by polymerization to many thousands of 
subunits [MK], The resultant MT polymer serves as a semi-rigid structural element in the 
cell that is critical for intracellular distribution networks, chromosome segregation before 
cell division, and neuronal activity. MTs within cells are typically initiated from a complex 
of proteins that form a nucleation template. Polymerization proceeds at a rate that is depen- 
dent upon the concentration of free tubulin subunits. Over time, a remarkable phenomenon 
is observed both in cellular MTs and with purified tubulin. The MTs undergo stochastic 
switching between states of polymerization and depolymerization, a process termed 'dy- 
namic instability' S(~!C'Q~TTj . Switching frequencies are modulated during certain cellular 
transitions to alter the length distribution and density of the polymer array. The mech- 
anisms governing dynamic instability are still an active subject of both experimental and 
theoretical investigation [KM]. 

In general terms, the dynamic instability arises due to an energy dependent mechano- 
chemical transition that occurs in the beta-tubulin half of the tubulin dimer after polymer- 
ization. MT polymerization proceeds in a head to tail fashion from the nucleation template 
leaving the beta-tubulin face exposed at the microtubule end. Free tubulin dimers rapidly 
bind guanosine triphosphate (GTP), a small molecule used as a convertible form of stored 
chemical energy in the cell. If the exposed beta-tubulin on the MT end is GTP-bound, 
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it has a relatively high affinity for binding the alpha-tubulin face of any free subunits in 
the cytoplasm. Once bound, the beta-tubulin undergoes a structural change increasing the 
probability of hydrolyzing the associated GTP to guanosine diphosphate (GDP). The bind- 
ing affinity of GDP-bound subunits for other tubulin dimers is relatively low. Therefore, 
GTP hydrolysis serves as a switch for changing the binding affinity of the subunit. If the 
rate of GTP hydrolysis in the MT outpaces binding of new GTP-bound dimers, the MT 
will lose its theorized cap of GTP-bound dimers and transit from a state of high affinity 
dimer binding to very low affinity. Over an important range of free tubulin concentrations, 
this loss of the GTP-cap will result in a switch from polymerization to depolymerization, 
termed 'catastrophe.' Should a GTP-bound dimer manage to bind the depolymerizing end, 
the polymerization can be recovered, a process termed 'rescue.' The observed stochastic 
switching between states of growth and shortening is attributed to the statistical variance 
in the processes associated with dimer binding and GTP hydrolysis. 

A closed form kinetic description of microtubule dynamics has not been successfully 
rendered to date. Several analytical models have been developed to describe how the major 
factors leading to dynamic instability (i.e. growth and shortening velocities together with 
catastrophe and rescue frequencies) will produce a steady state system of polymers under 
various conditions. The stochastic nature of the switching between polymerization states 
complicates the construction of determinative equations that can be numerically solved for 
properties such as MT length distribution. In this work, we advance the analysis of MT 
dynamics by extending the analytical methods applied to the most prominent model for 
dynamics instability in an in vitro context |DL| . 

2 Preliminaries 

2.1 Dichotomous Markov Noise 

The dichotomous Markov noise (DMN) v(t), is a simple two- valued stochastic process such 
that the state space of the random variable consists only two values {if 1 * 1 } with con- 
stant transition rates and /7. This means that the waiting-times in the two states are 
exponentially-distributed stochastic variables. The switches of v(t) are Poisson process with 
probabilities p~ and p + [IB]. We can describe these probabilities by a first-order kinetic 
equation: 
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Figure 1: Dichotomous Markov Noise. 



where p + (t) + p (t) = 1 for all times, and r c 
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-T= — -— — —r shows the mean time between 

^ f+ +/- 

switches of v(t) or the velocity relaxation time. The instantaneous average velocity is defined 
as: 



(2.4) 
Where 
(2.5) 
and 

(2.6) 



v(t) = v + p + (t) - v~p~(i) = v(oo) + (v(0) - v(oo)^je 



-Tt 



v(Q) = v + p + (0) -?Tp"(0) 



V oo 



This noise is a time-homogeneous Markov process and is therefore completely characterized 
by the following transition probabilty: 



P i3 {t) = Pr{v{t) = i\v{0) = j), i,j E {±^} 

The temporal evolution of the noise is given by the Kolmogorov forward equation, in the 
physical literature known as master equation [IB]: 



(2.7) 



dt \P +j (t) 



-n 



ft 



where k+ and k- are the mean frequencies of passage from A+ to —A- This system can be 
described by its transition matrix: 



(2i 



P__(t) P- + (t)\ _ f++fle 
P+_(t) P ++ (t)J Tc U+(i 



+ e- 
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In the stationary case we have: 

(2.9) Pr(v = v+) = f+r c Pr(v = -v-) = f+r c 

We also assume that the external noise is a stationary random process and hence the 
DMN has to be started with (2.9) as initial condition. The corresponding stationary mean- 
value is : 



(2.10) 



< v ( t ) >=V = (v+f+ - V~f+)T C 



2.2 Master Equations 

Neglecting the free tubulin concentration variations in the microtubule dynamics process, 
the dynamic instability equations governing the time evolution of density functions p + (x, t) 
and p~(x, t) are: 

(2.11) = _ f - p+{Xti) + fy- {x , t) 

(2.12) d ^=v- d J^ +r+ P + (*.t)-fy-( X , t ), 

where p (x, t) is the density function of growing (shrinking) microtubules in the interval 
[x, x + 8x\. (Hence, p jp is the probability density function of growing (shrinking) in the 
interval [x, x + 5x], where p := p + + p~ denotes total density.) 

In prior work of Dogterom and Leibler, the main result was the prediction of a sharp 
transition between an unlimited growth, with the average speed J > and a steady-state or 
bounded growth, characterized by a MT length distribution with J = 0. In the unbounded 
growth region the average length increase is < L >= Jt, where 

(2.13) j__v + r + -v-n 



and the distribution approaches asymptotically to a Gaussian of width s^yD e fft, where 

( 2 - 14 ) D *ff= i/i u f -^ v++v " 2 



and 

In the steady state the distribution of MT length is exponential with mean: 



v t) T 

(2.15) < L >= — + — 

v fl - v+f. 
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3 Model 



Following the Dogterom and Leibler model for growing/shrinking MT, we consider a general- 
ized one-dimensional concentration-dependent model for microtubule growth and shrinkage: 

d(v+p + (x,t)) . , , _, v d 2 (p + (x,t)) 

- ax - f + p + (^t) + f-p (m) + d yp d y 

., .. </// \-i.t) dp~(x,t) , . , i . . 1 d 2 (p~(x,t)) 

' ' ' = W + UP (X ' * } " ^ P ~ (X ' t} + ftc* 

- -fee x,t )+v-p- (x,t )-v + p + (x,t ) + D 



dp + (x 


t) 


dt 




dp~ (x 


t) 


dt 




dc(x 


t) 




Figure 2: Model geometry for infinitely rigid MTs originating from a planar surface 

In these equations, c(x,t) represents the concentration of the free tubulin, k is the 
hypothetic nucleation rate and 

(3.2) /_ = Loc(x,t) and v + = u + c(x,t) 

are the other concentration-dependent parameters. 

The constant k is assumed to be proportional to the number of available nucleation 
sites per volume; kc, measuring the rate of nucleation (with associated loss of free subunits) 
is thus proportional to the cross-section of a free subunit meeting a nucleation site, the 



simplest possible assumption. Likewise, (3.2) are derived from probabalistic cross-sections. 



Note that, properly, each of the coefficients u, u , and k should be assumed proportional 
to diffusion constant D measuring mean free path of subunits. As we hold D fixed in this 
analysis, there is no harm in taking them as constants however. 
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In this setup the nucleation points assumed hypothetically to be uniformly distributed 
throughout [0, +00). This model in the matrix form is as follows: 



(3.3) U t = A{U)U X + B{U)U + CU X 
on x > 0, where 

(3.4) U(x,t)= lp-(x,t) 

\ c(x,t) 

-u + c(x,t) — u + p + (x,ty 

(3.5) A(U) = I v- 



— /7_ ujc(x,t) 

(3.6) B(U) = I /" -uc(x,t) 



-u + c v —kj 



and 



(d 

(3.7) C = d 

\0 D 

and it + , fe, /+, w, d, D are positive constants. Here, ^ > are densities of growth/decay 
of microtubules and c > is concentration of free tubulin. 
Typical value for the parameters are 

D = 0.5, u = 0.15, v~ = 0.05, /" = 0.0005, d ~ 0, 

with u + , A; positive constants for which we don't know yet the typical values. 
Boundary conditions at x = are Dirichlet conditions on the densities p 1 * 1 , 

and on the concentration c either a Dirichlet condition 

(3.9) c\ x=0 = c 
or a homogeneous Neumann condition 

(3.10) d x c\ x=0 = 0. 
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3.1 Goals 



We seek to study the existence and stability of steady-state solutions of ( 3.3 H 3. 10) of 



"boundary-layer" type, that is, for which the solution U approaches a constant state U+ as 
x — > +oo. Precisely, we seek asymptotically constant stationary solutions 



(3.11) 



U(x, t) = U{x), lim U{z) = U+, 



of (3.3), i.e., asymptotically constant solutions of the steady-state ODE 
(3.12) A{U)U X + B{U)U + CU XX = 



with boundary conditions (3.8) and (3.9) or ( |3.10 ) at x = 0. When such solutions exist, we 
seek to study their stability, that is, whether a perturbation U satisfying the same equation 



(3.3) that is close to U at initial time t = in some norm remains close to U for all t > 



in some (possibly different) norm. In applications, it is only stable steady states that are 
truly steady in a practical sense; unstable steady states, though steady in an idealized sense, 
persist indefinitely only for a measure (hence probability) zero set of initial data. 

4 Calculation of the endstates [/+ 



Preparatory to the study of steady state solutions U of (3.11), we first investigate the 



possible limiting states U+ to which U may converge, seeking rest points 
(4.1) B(U+)U+ = 



of (3.12) 



The first two equations of (4.1) yield p 
(4.2) 



n 



U ' UJC 



+ u p 



■p + , whereupon the third equation yields 



ck. 



Solving, we obtain the one-parameter family of solutions 

kcf+ + ukc 2 



(4.3) 



e ■.= { P - 



P 



c > arbitrary} . 



V f , — U+UJC 2 V / 1 — ii+WC' 

From the requirement < p~ < +oo, we obtain the physicality condition 
(4.4) u~fl -u + ujc 2 > 0. 



Remark 4.1. Equations (4.2), (4.4) simply reflect the pair of balance laws f+p = f_p~ 
and v~p~ — u + p + = ck > expressing conservation of densities and free subunit concen- 
tration, respectively, in the absence of spatial variation. Recall that ck is the rate at which 
free subunit concentration is lost to nucleation, or binding of subunits to nucleation sites. 
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4.1 Exponential decay to Uj, 



We next investigate the rate of decay of solutions approaching endstates U+, specifically, 
whether or not approach at uniform exponential rate. As there is a continuous curve £ of 
viable endstates, this amounts to verifying that the linearized ODE about L7+ , written as a 
first-order system, has a center subspace of dimension one, with no other neutral modes. 



Linearizing (3.12) about a rest point U = U+, we obtain 
(4.5) 



B(U+)U + A(U+)U X + CU XX , 



where 



(4.6) B(U, 



'0 

B(U+) + | 
,0 



U = e» x V in (4.5), to obtain 




-f+ 
f+ 



UJC+ 




It is convenient to find eigenvalues of (4.5) directly in second-order form, substituting 



or 



(4.7) 



(pA(U+) + B(U+) + i?C)V = 0, 

= det(/M(£/+) + B(U+) + \i 2 C) 

(d/j, 2 — u + c + [i — uc + 
/7J d/i 2 + v~[i — uc+ 

— u + c + v~ 

(dfj, — u + c + dfj, + v~ 

/7J dpi 2 + ufi — wc-| 
— u + c + v~ 



— uip + — u^p+iJ, 
up^ 

c[i 2 — u + p\ — kj 

-u + p+ 
up+ 

- u + p\_ — 



(4.8) 

q{p) := (Dd 2 )fi 5 + (Ddv- - dDu+cj^ 

- [d 2 u + p\ + d 2 k + f+dD + u + c+u~D + duc+D^j 

+ ( — f^v~ D + u r c?,iijD + u + c + dk — dv~u + p\ — dv~k)fi 2 



+ ydujc+k — dup + v + u + c + u k + dojc+u^ p~\_ + f + dk + f + du + p^ 
+ (f^u^k — u + (? + ujk ) , 



u c + dcop , )/i 



where we have obtained the third equality in (4.7) by adding row two to row one and 
dividing out \i from the result. 

So long as q has no pure imaginary roots, then the center subspace is dimension one 
and we obtain exponential decay. We now check this, first for zero roots, then for nonzero 
imaginary roots. 
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4.2 Checking zero roots 



Plugging in \x = 0, we get /, 



-4- 2 



0, which is excluded by (4.4). Thus, we may 



conclude that there are no additional zero roots. 



4.3 Checking pure imaginary roots 



Plugging in // = i£ in (4.7) we get 
(4.9) 



(Dd 2 )if + \ Ddv~ - dDu + cjt + [d 2 u + P X + d 2 k + f+dD + u + c + u~D + dwc+DUf 
- ( - f±v~D + uc\ujD + u + c + dk - dv~u + p\ - dv~k\ 2 

+ (dujc+k — dujp+v~ + u + c+v~k + dwc-fU + p^ + f+dk + /^Td , u + p+ — u + c+(fwp^i£ 
+ ( f7v~k — u^c+uik 



it^Dd 2 ^ 4 + (d 2 u + p\ + d 2 k + f+dD + u + c + v-D + duc+D^ 2 
+ (duc + k — dujp+v~ + u + c + v~k + dujc + u + p\ + f+dk + f+du + p\ — u + c + dup^j^ + 
(^Ddu- - dDu + c^ A - ( - ftv~D + uc^W-D + u+c+d/c - - dv'kjf 



+ (/, z/ k — u + (? + ujk 



We need only check nonexistence of pure imaginary roots fi = i^ with £ 7^ 0, £ real. 
Noting that q(i£) = ^i(£ 2 ) + (/2(£ 2 )> with q\ and (72 quadratic, we find that such £ must be 
of form £ = ±y/x, where x is a common real positive root of q\ and 92- The unique common 
root is easily found by taking the resultant (taking a multiple eliminating the quadratic 
term and solving the resulting linear equation, then substituting this result into q\ or q^ to 
check whether it vanishes- also checking whether it is positive) . This computation is carried 
out in Appendix |Aj 

Note that the resultant procedure described gives in the end an analytic function of the 
parameters of the problem, which vanishes if and only if q\ and 52 have a common root. By 
properties of analytic functions, it therefore vanishes either on a surface of measure zero in 
parameter space, or else for arbitrary choice of parameters. But, it is readily checked for 
specific parameters that there is no common root (see Appendix [A]). Collecting information 
we may conclude the following general fact. 

Proposition 4.2. For generic choice of endstates U+ € £ (i.e., all but a set of zero one- 
dimensional measure), all orbits converging to U+ do so at uniform exponential rate. 
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5 Dimension of the stable manifold at U + 

We next determine the dimension of the stable manifold at [/+, that is, the number of 



eigenvalues [i of (4.7) with negative real part, or, equivalently, the number of roots of 



q(n) = 0. As a first step, computing the mod-two stability index 

sgnq(0)q(+oo) = sgn (f7u~k — u + c 2 + u:k){Dd 2 ) = +1 



using (4.4), we find that the number of unstable (i.e., positive real part) roots is even, so 
that the number of stable (negative real part) roots is odd, and thus is 1, 3, or 5. 

Next, using a standard technique in the stability theory, we study the dispersion relation 
for the time-evolutionary problem, from which we may conclude by homotopy argument that 
the dimension is three. Specifically, we look at the linearized time-evolutionary equations 
about U+, 

(5.1) \U + B(U+)U + A(U + )U X + CU XX , 

written as a first-order system, and look again at the number of negative real part eigen- 
values, or solutions ^ of the indicial equation 

(5.2) = det(A£7 + fiA(U+) + B(U+) + n 2 C) 

on the domain A := {A : KA > 0} of interest for the later stability analysis. 
Setting fx = i£, £ real, we obtain the dispersion equation 

(5.3) = det{-\U + i£A(U + ) + B{U + ) - fC), 
determining a family of three curves 

(5.4) X,(0 e a(iSA(U+) + B(U+) - £ 2 C), 

j = 1, . . . , 3, where cr(M) denotes the spectrum (eigenvalues) of a matrix M. 



We make the following standard assumption, verified in Section 5.1 for c+ small and 
checked numerically for the specific cases treated in this paper. 

Assumption 5.1. The constant solution U = {/+ is linearly stable, i.e., 
(5.5) &Aj(C) < 

for all £ G 3?, with equality only for £ = 0. 



Lemma 5.2. Under Assumption 5.1, the numbers of stable and unstable roots of fi of (5.2) 



are independent of A for all KA > 0, A ^ 0, and are both equal to three. 
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Proof. The roots are continuous as functions of A, so for the first assertion we need only 
check that no roots can cross the imaginary axis, i.e., there are no roots fi = i£ with £ 
real and 5RA > and A / 0. But, this possibility is ruled out by Assumption |5.1[ Thus, 
the numbers of stable and unstable roots are constant. To determine their values, we may 
take A —> +00 along the real axis, and compute directly that in the limit there are three of 
each. □ 



Corollary 5.3. Under Assumption 5.1, the dimension of the stable manifold at states 
U+ G £ is generically three. 



Proof. By continuity of the roots of (5.2) with respect to A, together with the fact already 
established that at A = there is generically a single zero root, we find, taking the limit as 
A —7- 0, that there must be for A = either two or three stable roots, depending whether the 
single zero root is the limit of a stable or of an unstable root as A — > 0. Since the number 
is odd, by our index computation, we conclude that the dimension is generically three. □ 

Remark 5.4. Note by our computations that we have obtained also the additional in- 
formation, which will be important for later stability analysis, that the zero root at A = 
perturbs as the real part of A is increased to the unstable half plane. As we shall see 
later, this has the important consequence that the single undamped mode governing "total 
density" p := p + + p- is convected inward toward the boundary x = 0, and not away 
toward x = +00, and therefore the nonlinear stability theory may be treated by simpler 
weighted-norm methods [He, SatJ rather than the delicate pointwise methods of |YZ1 INZ] . 

5.1 Verification for c + small 



We now verify Assumption 5.1 for c+ small, by explicit computation of the case c+ = 0. 



Proposition 5.5. Assumption 5.1 holds for U+ G £ and c+ sufficiently small. 



Proof. Plugging U = e^ x in ( |5.1[ ) and the fact that A G a(B + Ai£ - C£ 2 ), we get 
(5.6) 




-1i + C4 


—u + p\ 







v- 


) 












— uip+ — i^u 




ed 


ujp+ 






—k — u + p\_ — 


e 




For case £, plugging 

(5.7) p+= " kC \ , , * C/ + + 2 
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in ( |5.6[ ), we get 
(5.8) 

f -n - & + c+i - ed 

D+M<,-(t= f+ 

—u + c + 

When c+ = 0, the eigenvalues are 



V 



UJC-\ 



-ojc + + £v i — ^ 2 d 



ktoc-\-{f, +£u~T~a+i) 

— f — /T+if+tJC 2 , 
k{-v-f--t-Du-f-+ZDu+plc i + ) 



(5.9) 



A i 



-f7-ed, \ 2 = -S 2 d + i£v-, A 3 



which indeed have nonpositive real part for all real £. 

In the general straightforward computation shows that at £ = there is a single 



zero eigenvalue (note that the lower lefthand 2x2 block of B is nonsingular, by (4.4)) of B 
with left and right eigenvectors L = (1, 1, 0) and R = (a, b, 1) T , where a + b 



f+k 



v f + -uc%u 



2 „.+ 



and 
(5.10) 



a := LA+R 



2 -I 



a + b 



>0, 



whence, by standard matrix perturbation theory, the Taylor expansion of this eigenvalue at 
£ = is 

HO = iat + Pt 2 + '- + 5S 3 , 

a > verifying directly our observation of inward convection in the neutral mode. (Note: 
the fact that a is real follows with no computation, simply from ( 5.10| ).) Moreover, we may 
verify 3ff Ai (0) < 0, 9iA3(0) < by computing the characteristic polynomial 

det(l? + — fj,I) = /x(/i 2 + r/i + s) 

and verifying the discriminant condition 

r 2 - As = (u + p\ + fc + wc+-/;) 2 + 4(f~k + u + p+) > 0. 

For c+ sufficiently small, < 0, 8 uniformly bounded, by continuity of the Taylor 
coefficients in the limit. This establishes that KA2 < for £ near the origin, with equality 
only at £ = 0. On the other hand, 9iAi and 5RA2 by continuity are strictly negative on any 
bounded set of £ and all three are strictly negative (almost by inspection) in the limit as 
|£| — > 00. This establishes the assumption for c+ sufficiently small. □ 

5.2 Verification for general c + 



For general c+, we verify Assumption 5.1 numerically in the course of our Evans function 



computations, via the following elementary observation. 
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Lemma 5.6. For D,d > 0, Assumption \5.1\ is equivalent to the property that for A pure 
imaginary and < |A| < R, R > sufficiently large, the number of roots \jl of indicial 



equation 5.2 with negative (resp. positive) real parts is constant and equal to three. 



Proof. It is evident that for |£| sufficiently large 9iAj(£) < — #£ 2 , some 9 > 0. Thus, if 



Assumption 5.1 is violated, then 3iAj(£) is pure imaginary for some < |£| < R, or, 
equivalently, there is a pure imaginary root [i = i£ of the indicial equation |5.2| for some pure 
imaginary A with < |A| < R. But, this means that either the number of negative real part 
roots or the number of positive real part roots must be less than three, or else the total of 
all roots would be at least seven, a contradiction. □ 

As described below, the numbers of stable and unstable roots [i of the indicial equation 
are checked at each step of our Evans function computations, in particular, on a n im aginary 



interval [—iR, iR] with R larger than the value R described above; see Remark 8.2 



6 Existence of steady-state solutions 
6.1 General theory 



Writing (3.12) in first-order form as 

(fU) \p) x = - B(U)U 

with P := U x , we obtain a six-dimensional dynamical system. The three boundary condi- 
tions at x = determine a four-dimensional manifold of solutions (three free parameters 
in the initial value, plus one dimension in the direction of spatial-evolution for a particular 
initial value at x = 0). 



By Corollary |5.3[ under Assumption 5.1 the stable manifold associated with a rest point 
[/+ generically has dimension three. Thus, in looking for stationary solutions satisfying the 
boundary conditions at x = and converging as x — > +oo to a specified endstate U+, 
we seek the intersection of a four-dimensional with a three-dimensional manifold, which 
should generically consist of a one- dimensional manifold, or a finite union of distinct curves, 
corresponding to distinct steady-state solutions. 

That is, by a dimensional count, connections seem to be possible for all endstates U+ € £ 
that are stable as constant solutions, as (numerically) all examples considered seem to be. 

Remark 6.1. The above suggests, by formal matched asymptotic expansion, that behavior 



far from the boundary (distance >> max{d, D}) of solutions of (3.3) should be governed 
by the inviscid equation Ut = A(U)U X + B(U)U, with boundary condition U G £ at x = 0, 
whether for Dirichlet or Neumann conditions imposed in the full viscous equations, so long 
as the boundary layer (i.e., steady state solution connecting to £) is stable. See [GMWZ 
for related discussion in more general context. 
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6.2 Numerical determination of steady-state profiles 



To determine the profiles U of steady-state solutions, we use a numerical boundary- value 
solver, imposing boundary conditions at x = and projective boundary conditions at +00, 
ensuring correct entry toward a specified endstate in £ along its stable manifold. 



Expanding (3.12) we have, 



(6.2) 



= —u + cp+ — u + p + c x — f + p + + cocp + dp£ x , 
= v~Px + f+P + ~ ucp~ + dp~ x , 
= — u + cp + + v~p~ — kc + Dc xx , 

which yields the first order system 



(6.3) 

with 
(6.4) 



fyi\ 
2/2 
2/3 
2/4 

2/5 

W 



F(Y) :- 



/ 2/2 \ 

l{u + y 5 y 2 + u + yiy 6 + /+2/1 - ujy 5 y 3 
2/4 

3(^2/52/3 -v~y±- f+yi) 

2/6 

^(« + 2/52/i - ^~2/3 + ky 5 ) 



^ = (2/1,2/2,2/3,2/4,2/5,2/6) -=(p + ,pt,P ,Px> c > c xY 



The Jacobian at +00 is given by 



(6.5) 



dF 
W 



( 

d 


d 


V-' 



1 

u + c + 
d 








D 





WC-l- 

~~d~ 






D 






1 

_ v2_ 
d 








d 





Z5 



\ 

d 




1 

y 



withy + = (p+,o,p;,o, c+ ,o) T 

6.3 Numerical method 



We solve (6.3) as a two-point boundary- value problem on [0, M] with M > chosen suffi- 
ciently large. At x = 0, we impose the three boundary conditions 

(6.6) yi(0) = pj, 2/3(0) = Pq , and 2/5(0) = c or y 6 (0) = 0. 

At x = M, following the general approach of |Bej . we impose projective boundary conditions 

(6.7) L ■ ( yi (M) - pX) = L ■ (y 3 (M) - p~) = L • (2/5 (M) - c+) = 0, 
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where L is any 3x6 complex matrix that is full rank on the center unstable subspace of the 
Jacobian ^y(Y + ). We use the simple and well-conditioned choice of L consisting of rows 
spanning the unstable left eigenspace of |^(Y + ). 

Following [BHRZt lHLZjlCHNZj . we implement this scheme using MATLAB's boundary- 
value solvers bvp4c, bvp5c, or bvp6c [HMJ, which are adaptive Lobatto quadrature schemes 
and can be interchanged for our purposes. The value of approximate spatial infinity M are 
determined experimentally by the requirement that the absolute error \U(M) — U+\ be 
within a prescribed tolerance, say TOL = 10~ 3 . For rigorous error /convergence bounds for 
these algorithms, see, e.g., |Bej . 



6.4 Results 

For each choice of endstates and parameters tested, we successfully found steady-state 
profiles. A typical example is shown in Figure 3 for endstate Z7+ = (0, 0, 0) and parameters 
d = .1, D = 1, U3 = 1, u~ = 1, /7 = 1, u + = 1, k = 1. Other examples may be found in 
Figures 7-9, superimposed on time-evolution plots for nearby perturbed solutions. 

7 Abstract stability theory 

Under the assumption D,d > 0, C is positive definite, so (3.3) is a system of semilinear 
second-order parabolic equations of the type considered in |He} ISatj . Stability of steady- 
state solutions of this system may be treated in straightforward fashion by the weighted 
norm method of [SatJ, as we now briefly describe, reducing the problem of determining 
nonlinear stability to that of checking spectral stability of the linearized operator about the 
wave. The simplicity of our treatment here is to some extent an accident of the favorable 
structure of this specific problem. For (more complicated) methods applying to general 



systems, see, e.g., |HZ1 IZ2j and references therein. 
7.1 Abstract linearized equations 



Linearizing (3.3) about a steady state U(x), we obtain the linearized equations 

(7.1) U t = LU := A(x)U x + B(x)U + CU XX , 

A{x) := A(U(x)), B(x)V := B(U(x))V + (dA(U)V)U x , with boundary conditions 

(7.2) (p + ,p-,c) = (0,0,0) or (p + ,p-,c x ) = (0,0,0). 

7.2 Abstract eigenvalue equations and spectrum 



The eigenvalue equations associated with (7.1) are 



(7.3) (L - X)u = 
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Figure 3: Sample profile: (p+,p + , c+) = (0,0,0), d 
u + = 1, k = 1. 
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with boundary conditions (7.2). Let W '°° as usual denote the Banch space of bounded 



measurable functions on x G [0, +oo) possessing a bounded measurable weak derivative, 
with norm [|/||^i,oo := sup|/| + sup|/'|. Following |Hej . define the spectrum a{M) of an 
operator M with respect to a given Banach space B to be the set of A € C for which (A — M) 
does not possess a bounded inverse on B. Define the point spectrum a p {M) of M to be the 
set of A G C for which (A — M)U = has a nonzero solution in B, and the essential spectrum 
o~ ess (M) of M to be the set of A G C that are in a{M) but not in a p (M). 

Lemma 7.1. Under Assumption 5.1, 5?A < for any A G a ess {L), where spectrum is 
defined with respect to W ,ao , with equality at A = 0. 



Proof. By Lemma 4.2, the coefficients of L converge asymptotically as x — > +oo, whence, by 
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a standard result of Henry |He| (adapted in straightforward fashion to the case of the half- 
line), the essential spectrum of L is either bounded to the left of the rightmost envelope 



of the dispersion curves Aj(£) of (5.3) or else contain all points to the right. The latter 



possibility is easily eliminated by a standard energy estimate, whence the result follows by 



(5.5) and the fact that A2(0) = (computed earlier). □ 



7.3 The method of weighted norms 

Following [Sat], introduce now the weighted norm 
(7.4) [|17||„ := \\e^ x U\\ w x,oo, 

rj > 0, and the associated Banach space B of functions with bounded || • [L norm. The fol- 
lowing key lemma states that this weighting has the effect of shifting the essential spectrum 
of the linearized operator into the strictly stable half-plane {A : KA < 0}. 



Lemma 7.2. Under Assumption 5.1, for rj > sufficiently small, KA < —0{r}), 9 > 0, for 



A E cr ess (.L), defined with respect to Banach space B . 

Proof. Making the change of variables U — > e vx U following [Sat], we convert L to the 
operator 

(7.5) L := C(d x - r,) 2 + A(d x - rj) + B =: Cd 2 + Ad x + B, 

where B = B + r]A, etc. (we need not compute A and C). Checking definitions, we find 
that the spectrum of L with respect to B is the spectrum of L wit h re spect to W 1,co . 



Applying the standard theory of |Hej as in the proof of Lemma 7.1 we find that o~ ess (L) 
with respect to W 1,co is bounded on the left by the rightmost envelope of the dispersion 
curves 



Xj(0 G a(B + + i£A+ - eC+) = <j[B + + (ii- V )A + + (i£ - V ) 2 a 



Observing that Xj are continuous in rj, we find by Lemma 7.1 that, for 1/R < |£| < R, 
R > arbitrary, JRAj < —6, 6 > for n > sufficiently small. Likewise, < —9 for |£| 
sufficiently large and r\ sufficiently small, by domination of term C£ 2 . 

Referring now to the computations in the proof of Proposition |5.5[ and following the 
notation therein, we recall that JRAi and KA3 are strictly negative for all £, so that the above 
continuity argument in fact gives 3?Ai, A2 < —0 for r] > sufficiently small. 

Hence, it remains only to check the behavior of A2 for |£| < 1/R arbitrarily small. 



Recalling the Taylor expansion (5.10) of A2(£), and substituting A2(£) = A2(£ — rj/i), we 
find that 

KA 2 (0) = ma{i - rj/i) + h(£ - rj/i) 2j ) + 0(|£ - n/i\ 2j+1 ) = -ar) + K/i£ 2j + o(\£\ 2j + rj), 
where a := LAR > and h is the first term with nonzero real part in the Taylor expansion 



of A2 after a, by Assumption 5.1 necessarily an even-order term with K/i < 0. Thus, 



KA2 < —9 < for some 9 < 0, for \rj,^\ sufficiently small, completing the proof. □ 
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7.4 Basic nonlinear stability theorem 

We make a final assumption on stability of point spectrum of L. 
Assumption 7.3. 9?<t p (L) < 0. 
Theorem 7.4. Under Assumptions 



5.1 



and 



7.3, for rj > sufficiently small, let U(x,t) be 
a solution of (3.3) with initial data Uq such that He^tZo — J7)||jyi,oo is sufficiently small. 
Then, \\e vx (U — U)\\wi,<x>(t) decays to zero as t — > oo ; at rate 



(U-U)\\ w i,oo(t) < Ce 



-et 



where 9 = 9{rf) > 0. (Here, 9(rj) ~ an as rj — > 0, where a is as in (5.10).) 



Proof. By Assumptions 5.1 and Lemma 7.2 we have Kcr ess (L) < —9 < with respect to 
space B. On the other hand, eigenvalues of L with respect to B are necessarily eigenvalues 
with respect to W 1,co as well, hence, by Assumption 7.3 Ko" p (L) < —9 < with respect to 
B as well, and so 



(7.6) 

with respect to the space B. 



sfto-(L) < -9 < 



By the spectral gap (7.6) and the fact that L as as second-order elliptic operator is 
sectorial, L by standard analytic semigroup theory [He] is linearly exponentially stable with 
respect to £>, i.e., 



f\\v < e~ 



Turning to the nonlinear equations, and noting that e vx is an exponentially growing weight, 
we find as in [SatJ that nonlinear exponential stability follows as well, by a standard nonlin- 
ear iteration/ variation of constants argument. We omit these straightforward details. □ 

Theorem |7.4| asserts that, similarly as for finite-dimensional ODE, exponential stability 
of steady-state solutions of the PDE (3.3) reduces to spectral properties of the linearized 
operator L: specifically, verification of Assumptions 5.1 and 7.3 This abstract result reduces 
the study of PDE stability to the study of the eigenvalue equation for L, an ODE. 

Remark 7.5. The success of the weighted norm method relies on the property that the 
undamped mode p = p + + p~ is convected inward toward the boundary x = 0. Intuitively, 
we may think of p to lowest order as p(x + at) for a fixed profile p, so that 



sup e vx p 



e~^ at supe 



decays time-exponentially for rj > 0. More general situations may be treated by pointwise 
methods as in |HZ|, IZ2| . to obtain time-algebraic rather than exponential rates of decay. 
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8 Spectral stability analysis 



Theorem |7.4| reduces the study of nonlinear stability to verification of spectral stability 



assumptions 5.1 and 7.3 the first a somewhat standard linear algebraic problem and the 
second a type of ordinary differential boundary- value problem arising frequently in the study 
of eigenvalues of differential operators. We treat these numerically using numerical Evans 
function techniques developed in \Bi\ IBDGl IHuZ| . 



8.1 Linearized eigenvalue equations 



Expanding (7.1) coordinate- wise, we obtain the linearized equations 

dtp + = —u + cp x ~ — u + c x p + — u + cpx — u + c x p + — f+p + + ucp~ + ujcp~ + dp xx 
d t p' = v~p~ + f~p + - ucp- - Lucp~ + dp xx 
dtc = —kc + v~p~ — u + cp + — u + cp + + Dc xx , 



1.1) 



with boundary conditions (p + ,p~,c) = (0,0,0) (Dirichlet) or (p + ,p~,c x ) = (0,0,0) (Neu- 
mann) . 

Seeking normal modes U(x,t) = e A *u(x), we obtain the linearized eigenvalue equations 



-u cp x 



1.2) 



U'CxP' 



Xp 
Xp 

Ac = — kc + v~p 
with boundary conditions 

(8.3) (p + ,P~,c) = (0,0,0) or (p + ,p-,c x ) = (0,0,0). 



V Px + f+P' - UCp 



u cp 



u cp x — u c x p 
' - ujcp- + dp xx 
' — u + cp + + Dc x 



f+ P + + WC P + ^cp + dp x 



We seek growing or neutral modes, i.e., bounded solutions of (8.2)-(8.3) with 5RA > 



8.2 The Evans function 



Written as a first-order system, the eigenvalue equations (8.2) appear as 

W x = A{x,X)W, 



(8.4) 
where 
(8.5) 
and 



1.6) 



W = (w 1 ,w 2 ,w 3 ,w 4 ,w 5 ,w 6 ) T := {p + ,pt,P ,p x ,c,c x ) T 



A 
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d d 
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x G [0, +00), with boundary conditions (8.3) imposed at x = 0. 
The limiting coefficient matrix as x — > +00 is 



(8.7) 



A+(A) := 



/ 

a+/; 

d 


d 



\ D~ 



1 

— d~ 













A+UIC4 

dT~ 








d 





D 








u+p\ 




d 














1 







) 



which for A = reduces to the Jacobian ||I(Y + ) computed in (6.5). By exponential con- 
vergence of U(x) to U + as x — > 00, we have exponential convergence also of A(x, A) to 
As- 
under Assumption |5.1[ 



By Lemma 5.2 and Corollary 5.3, under Assumption |5.1[ the stable subspace of A 
has dimension three for 3ftA > 0, whence, by a standard lemmgQ using the exponential 
convergence of A to A + , there exist choices of three independent solutions Wj~ , j = 4, 5, 6 
of (8.4) spanning the set of solutions that are decaying as x — > +00 that are analytic as 
W 1,oo [0,+oo). Likewise, one may easily construct analytic choices of 
j 



1, 2, 3 of (8.4) satisfying the boundary conditions (7.2) 



functions from x 
three independent solutions W® 
at x = 0. Evidently, A is an eigenvalue, i.e., there exists a solution of (8.4) that is decaying 
as x — > +00 and satisfies boundary-conditions (7.2) at x 
dependency among W®, W®, W° , W 7 "/, W§ 
Defining the Evans function 



0, if and only if there is a linear 



W f 



6 ' 



E(X) :=det«,Ty 2 °,W 3 ,^ 



4 ; 



w£,wj 



following |AGJ1 IGZ[ IZ1| , we thus have that the eigenvalues of L with respect to the weighted 
space B on KA > correspond to the zeros of E. This can be efficiently computed numer- 
ically using the STABLAB package developed by J. Humpherys, based on algorithms of 
|BrZ| IHuZj . We refer the reader to [HuZj |HLyZ| for a discussion of theory and numerical 
protocol, which is by now standard^] 



8.3 High-frequency bound 

Define 5 := mm{d, D}, a = m&x\A(x)\, (3 = max |i?(x)|, A, B, C the coefficients of L, where 
\M\ denotes the Euclidean matrix operator norm of a matrix M, i.e., the square root of the 
largest eigenvalue of M T M. 

Lemma 8.1. There exist no spectra of L with respect to B for 3ftA + |9A| > a 2 /5 + (3 and 
n > sufficiently small. 

1 The "gap" or "conjugation" lemma; see, e.g., [GZl IZll ETMWZ] . 

2 See, e.g.. rBHRZIIliT^ICHNZI[lTr^lEii^lBl^llBLeZ| . 
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Proof. First note that eigenvalues of L with respect to B are eigenvalues with respect to L 2 
as well. Taking the real part of the L 2 complex inner product of U with eigenvalue equation 

(8.9) XU = AU X + BU + CU XX , 

we obtain after an integration by parts of term J U*CUdx the inequality 

UX\\Uf L2 < a\\U\\ 2 L2 + P\\U\\ L2 \\U X \\ L2 - 5\\U x f L2 . 

Likewise, taking the imaginary part yields |9A| ||C/||^ 2 < ct \ \ U\\ 2 L2 . Summing, and using 
Young's inquality, 2aab < (a 2 /5)a 2 + 5b 2 , yields 



(SftA+|QfA|)||C/||| 2 < (a 2 /5 + p)\\U\ 



L 2 ' 



yielding + \$sX\ < a 2 /5 + f3 whenever || L 2 / 0. 

This establishes nonexistence of point spectra for 3f?A + |9A| > a 2 /5 + (3. Nonexistence of 
essential spectra follows by the same argument applied to the linearized eigenvalue equations 
about the constant solution U = U+ together with the observation [He] that essential 
spectrum of L with respect to 1? is bounded on the left by the rightmost envelope of 
the spectrum of the linearized operator about U = U+, given by the solution set of the 



dispersion relation (5.3), and the fact that nonexistence of spectrum with respect to B is 
implied by nonexistence of L 2 spectrum provided r\ > is taken sufficiently small. (See the 



argument of Lemma 7.2 ) □ 



Remark 8.2. The above energy estimate implies in passing that R := a 2 /5 + f3 > R, 



where R is as in Lemma 5.6 To verify Assumption 5.1, therefore, it is sufficient to check 



that the stable subspace of A + has dimension three for A on imaginary interval [— iR, iR] . 
8.4 Numerical stability computations 

Following a standard strategy introduced by Evans and Feroe [EF] , we may now check by a 
single numerical winding number computation both of Assumptions |5.l| and |7.3[ Specific ally, 



.1 



we compute the image under E of a semicircle S of radius R := a 2 /5 + /3 (see Lemma 
for definitions) in the unstable complex half-plane KA > 0, with diameter [—iR, iR] lying 
along the imaginary axis. 

As the code for approximating E in the course of the computation checks that the stable 



subspace of A + has dimension three, we obtain by Remark |8.2| a check of Assumption 5.1 



in the course of computing E on the diameter A G [—iR,iR]. Assuming that Assumption 



5.1 is valid, we then obtain by Lemma 8.1 that there are no eigenvalues of E outside the 



semicircle. 



Moreover, by the properties of the Evans function described in Section 8.2 (recall: also 



dependent on Assumption 5.1), eigenvalues of L within the semicircle S are exactly the 
zeros of the analytic function E. By the Principle of the Argument, the number of zeros 
within S is equal to the winding number of E(S), i.e., the number of times E(S) circles the 
origin in counterclockwise direction as A traverses S in counterclockwise direction. Thus, 
Assumption\7.3\ corresponds to winding-number zero. 
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8.5 Results 



Figure 4 displays the image of the semicircle S under the (analytic) Evans function E for 
the typical profile U displayed in Figure 3. We see clearly that the winding number is zero, 
indicating spectral stability, hence, by the results of Section [7| U is nonlinearly stable. Here, 
d = .1, D = 1, oj = 1, v~ = 1, /JJT = 1, u + = 1, k = 1, so that a = I, (3 = 1.4343, and 
5 = .1, and -R = a 2 /<5 + /3 ~ 11.43; the radius of S is taken as i? = 12 > i£, following Lemma 



8.1 Though we did not carry out a systematic study over all parameter values and profiles, 



for all profiles checked, we obtained results of zero winding number consistent with stability. 

Remark 8.3. Stability of "small-amplitude" steady states U with max \ U— E/+ 1 sufficiently 
small, could in principle be determined by a study of stability of constant solutions, as 
described in a more general setting in [GMWZ]. However, this does not seem particularly 
interesting for applications, and so we do not attempt to carry out such an analysis here. 



9 Behavior for general initial data 

The results of Sections [7] and [8] indicate time-exponential stability of steady state solutions 
U with respect to spatially-exponentially decaying initial perturbations. To put this an- 
other way, solutions with initial data converging exponentially as x — > oo to U+ £ E and 
sufficiently close to the profile U, will converge as t — > oo to U. 

However, as we explore in this section, the actual stability properties appear to be 
considerably stronger. Namely, initial data converging as x — > oo to any limiting value Uoo 
appears to settle down rapidly to a steady state, determined by the total density p := p + +p~ 
of the limiting state Uoo- (Indeed, we suspect that only p need converge to a limiting value 
in order to determine the limiting steady state.) 

The weighted norm methods used above do not appear sufficiently fine to establish this 
property, or at least we were not able to carry out such an analysis. An interesting open 
problem would be to attempt to establish this result using the more detailed pointwise 
Green function methods of [HZ, Zl]. 



9.1 Numerical time-evolution approximations 

In this section, we describe the results of a numerical time-evolution study, that is, the 



evolution in time of given initial data under equation (3.1). This was carried out via 
MATLABs Finite Element Method (FEM) routines, with Neumann boundary conditions 
at a boundary x = 50 far from the computational domain of interest. For reference, see the 
helpful notes |H]. 

In Fig. 5, we have displayed graphs of variables p + , p~ , and c taken at successive time 
intervals, with boundary data p$ , p$ , and Co at x = fixed at the same values .2, .2, .2 
used in the computations used to obtain the steady-state solution from Fig. 3, and step 
initial data consisting of .2, .2, .2 for x G [0, 3] and 0, 0, elsewhere. We use the same 
model parameters as for Fig. 3 as well. 
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Figure 4: Typical winding number computation, corresponding the the profile of figure 3. 
the displayed contour is the image under D of a semicircle of radius 12, whereas R « 11.43. 

We see that the solution rapidly smooths and settles to an apparently steady state 
represented by the lower envelope of the graphs, with the darker bands near the lower 
envelopes indicating convergence of each variable toward these curves. In Fig. 6, we have 
superimposed onto Fig. 5 the steady-state profiles that were displayed in Fig. 3, computed 
as solutions of a two-point ordinary differential boundary problem, showing a nearly exact 
correspondence between the apparent limit of the time-evolution PDE solution and the 
theoretical steady state ODE solution computed by completely different methods. 

This is strong evidence for the accuracy of both computations, and also for stability 
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of the steady-state solution as computed by yet a third different method in our numerical 
Evans function computations. 




Figure 5: Time evolution of typical data. 
9.2 Selection of limiting profile 

The intuitive explanation for the behavior described in the previous subsection is that only 
the special mode p := p + + is "undamped", with other modes decaying exponentially 
toward their equilibrium values; see the discussion in Section [4] of stability of constant 
solutions. Thus, even if the limiting value [Too of the initial data is not on the equilibrium 
curve 8, it will rapidly adjust dynamically toward a value C/+inf£ with the same total 
density p. 

This is why the data from Fig. 5 converges toward a zero endstate profile. In Fig. 7, 
we show the results of a similar experiment with the same boundary data and the same 
initial data for p^, but c initially constant and equal to .02; as predicted by theoretical 
considerations, we see convergence to an identical zero endstate profile. In Fig. 8, we 
display the results for initial data c = .02 and p initially zero as x — > +oo, but p^ nonzero 
as x — > +oo (to do this, we assign p~ a nonphysical negative value, but mathematically this 
still makes sense). 

Further confirmation is given by Fig. 9, in which we display behavior for initial data 
converging to nonzero value of p; again, the agreement with prediction is exact. 
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[ 





Figure 6: Comparison with steady state profile. 

10 Conclusion 

We have applied new mathematical tools to the question of how microtubule dynamic 
properties contribute to microtubule array construction. Our results show the complex 
nature of the predicted end-states for a relatively simple model of microtubule dynamics 
including nucleation rate and tubulin concentration. We emphasize that the mathematical 
tools introduced to handle he evident stochastic processes here are of general application, 
and should find use in related situations. Originally developed to study stability of shock and 
boundary layers in gas dynamics, magnetohydrodynamics, and viscoelasticity (See [BHRZ, 
iHLZl ICHJNZl |HLyZj iBHZl iBLZl IBLeZj l. these methods are capable of handling systems of 



essentially arbitrary complexity. It is our hope that these methods and their application 
will open the way to the study of more complicated and realistic Biological models. We 
view the present study mainly as a feasibility study for the application of new mathematical 
tools transferring technology originating in the study of continuum mechanics to the area 
of biological modeling. 

An important further direction from both Mathematical and Biological point of view 
is the study of the singular perturbation limit as d — > 0. In our numerics, we have mostly 
investigated the nonphysical case d ~ D conducive to good numerical conditioning, whereas 
d « D in applications. An important direction for further investigation, discussed in 
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Figure 9: Selection of endstate (iii) nonzero limiting concentration. 



preliminary fashion in Appendix [Bj would be the rigorous analysis of the singular limit 
d — > 0. This should in principle be possible using the same body of techniques applied 
here; see [ZlJ for an analysis in a similar spirit. As discussed at the beginning of Section |9j 
another very interesting open problem would be to establish stability of steady state profiles 
with respect to perturbations decaying exponentially only in the total density p and not in 
other modes. This should be accessible by the techniques of [HZ> IZlj . 
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A Resultant computations 

A.l Resultant Procedure for c + = 

(A.l) 



qi (x) = [Dd 2 y 2 + (d 2 u + pX + d 2 k + f±dD + u + c + u'D + duc + D)x 
+ (duc + k — dujp^u~ + u + c + v~k + dujc + u + p\_ + /+dA; + f+du + p\_ — u + c + dup^j 
Q2{x) = ^Ddv~ — dDu + cjx 2 — ^ — f^v~ D + u + c^uD + u + c + dk — dv~u + p\_ — dv~k^x 
+ ^ — f+v~k — u + c\ujk^ 

For c + = these two reduce to 

(A.2) 9i (*) = (Dd 2 ^x 2 + (d 2 k + /±d£>)s + (/+dfc) 

g 2 ( x ) = (iJdz/^ar 2 + (f±v~D + dv~k s jx+ ( - f+v~kj 

Subtracting times the second equation from 1/d times the first yields a constant, 

2f~k > 

from which we may conclude that q\ and qi have no common root, in particular, q does 
not have any pure imaginary roots. 

A.2 Resultant Procedure for c + > 

Assuming D = d = 1, v~ = k = 1, = 0, u + = and wc + = wp+ = wp+ = 1, we get 



(A. 3) qi{x)=x 2 + 2x 

Q2(x) = x 2 + X 

which don't have a common positive root. 

B The singular perturbation limit 

In practice, it is important to take into account the singular perturbation structure imposed 
on the problem by the small parameter d. Otherwise, numerics will be destroyed when we 
try to go into the physical parameter range d « 1, for one thing. For another, there is an 
advantage to singular limits, in that they tend to reduce to problem to composite problems 
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that are sometimes explicitly solvable. Finally, d is often set to zero in applications, and 
it is important to justify that this approximation is valid. In this appendix, we initiate 
the discussion by a singular perturbation study of the profile existence problem. A parallel 
study of the stability problem should be carried out, but appears to be more difficult. 

First observe, if d = 0, then the number of boundary conditions should be only two at 
the boundary x = 0: one for c, which satisfies a parabolic equation, and one for the single 
incoming mode among hyperbolic . So, we have too many boundary conditions for the 
slow problem d = 0. 

We may conclude from this that there is a "true" (i.e., in general nonconstant) boundary 
layer of thickness ~ d matching full parabolic to hyperbolic-parabolic boundary conditions. 
This satisfies the "fast" problem obtained by rescaling x — > x := x/d, and involves second 
and first-order derivs. but not zero-order terms. That is, the p ± equations become in x 
coordinates, setting d = 0, just 

Vtx = (" + P + )x, Pxx = ~{v~P~)x, 

a pair of conservation laws, and constant-coefficient as well. Integrating up, we get 

vt = v+ {p + -j> + (+oo)), p~ = -u~(P~ p~(+oo)), 

from which we find p + = constant, but p_ (+oo) is arbitrary, with an exponentially decaying 
solution. The c equation degenerates on the other hand to c = constant. 

So, the correct boundary data for the slow problem are p + (0) prescribed (from original 
boundary conditions), c prescribed (from original boundary conditions, either Dirichlet or 
Neumann type preserved) . 

Then we solve the slow problem obtained by setting d = in the original coordinates: 

(B.i) o = _^v^n _ + f+p-( X ,t) 

(B.2) = v~ dj -^ + f+P + (x,t) - f±p-(x,t) 

(B.3) = -kc(x,t) + v-p-(x,t) - iy + p+(x,t)+D d2c } X ' t \ 

Observing (by addition of the first two equations) that 
(B.4) u + p + — v~p~ = a, 

a constant, we may rewrite these equations after some rearrangement as 

' ojc fl \ , , uca 



{^P^)x = — - — r- z^P + ) — 

\u u + cJ V 



(B.5) 

( a\ k ( a\ 
{ C+ k)xx=D{ C+ k)- 
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together with (B.4). Note that (B.4) already selects the value p (0) as a function ofp + (0) 



and c(0) ; in principle determining completely the solution with no further computation, 



provided one exists. The rest of the analysis consists in verifying that (B.4) is indeed 
consistent with existence. 

The second, constant-coefficient equation together with the requirement that c be bounded 
as x — > +00 may be exactly solved as 



c{x) + a/k 



-xJk/D 



(c(0) +a/k), 



a arbitrary, since only the exponentially decaying mode is allowable. In particular, we find 
that 



(B.6) 



c(+oo) 



From (14.41), we find that the coefficient 



= —a/k. 



of 



in the first equation is 



negative as x — > +00, so the stable manifold has a mode in direction p + . 

Note that, by this observation, the dimension of the phase space for the slow problem 
is four, while the dimension of the stable manifold about a rest state is two (one decaying 
direction for c, and one decaying direction for p + ) and the dimension of the manifold of 
solutions satisfying the (two) initial conditions is three (two free parameters plus direction 
of spatial evolution). Thus, as for the full problem, the dimension of the intersection is 
generically a union of finitely many curves/steady state solutions; in particular, we expect 
a connection to each possible endstate. 

We now investigate in detail the possible solutions. There are two cases. 

Case I. (c(+oo) = 0) In this case, also a = by (B.6), whence p~(+oo) = by (B.4) 

and ^ + (+oo) = u + c(+oo) = 0. Noting that in this case c(x) = e~ x ^ k ^ D c(0), we may 



rewrite the first equation of (B.5) as 




Px = 

to see that also p + — > as x — > +00, as the solution of a homogeneous first-order scalar 
equation with negative coefficient near +00. Moreover, (p + ,p~) go to zero superexponen- 

— Xyfk I D 

tially, as e e , while c goes to zero only exponentially. Such a slow solution exists for 

any initial data p + (0), c + (0), hence a corresponding matched boundary-layer solution exists 
for any Dirichlet data p + (0), p~(0), c + (0). 

Neumann data: in this case, the only bounded solution of the c equation satisfying the 
boundary conditions seems to be the trivial solution c = 0, from which (B.4) gives p~ = 0. 



But then (B.l ) yields p + = as well, i.e., only the trivial solution is possible, and we cannot 



generate a slow solution of this type for general initial data. 

Case II. (c(+oo) > 0) Again, only the decaying mode of (c + a/k) is possible, so that 
the limit as x — > +00 of c is c(+oo) = —a/k, or 



O: 



-kc(+oo) < 0. 
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From (B.4), we obtain then p (+00) = [cjv )(u + p + (+oo) + k), and the first equation of 



( B.5h yields at +00 the equilibrium condition p + = — „^ kc - — r. Solving for p_(+oo) using 

1 1 U~ f, —U+UJC 1 



(B.4), we obtain the full conditions for £ , as we should. Likewise, by similar computations 
as above, we obtain a solution for each initial data and endstate. 

(Note: this is not a "generic" result as in the general case. In other words, the assump- 
tion d — > precludes nongeneric behavior.) 

Neumann data: Again, the only bounded solution of the c equation satisfying the bound- 



ary conditions is the trivial solution c = —a/k, greatly simplifying the analysis. From (B.4) 



we obtain then p = (cjv )(u + p + + k), and the first equation of (B.5) becomes 



ujc 2 u + — f , V- \ , ojck 



Px I v~u + c ) P + V ll~ 
yielding at +00 the equilibrium condition p + = — " fcc , — 7 . Solving for p„(+oo) using 



(B.4), we obtain the full conditions for £, as we should. Assuming (4.4), we again obtain 
a connection for each endstate in £ , as expected. So, this case appears to permit profiles, 
also for Neumann data. 

Remarks B.l. 1. Neumann boundary conditions do not seem consistent in case a = 
corresponding to c+ = 0. Thus, this solution may (for Neumann conditions) be expected 
to disappear in the singular limit d — > 0. 

2. At a heuristic level, one could just drop one boundary condition and set d = 0, 
it appears, without much difference in profile behavior. However, for stability studies one 
should not neglect the fast (inner layer) dynamics near the boundary. 

In principle, stability should likewise be determinable from the singular perturbation 
structure in the limit as d — > 0, similarly as in |Zlj . However, we do not attempt to carry 
this out here. Some such analysis should be carried out in justification of the approximation 
d = commonly used in practice. 



35 



